home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / SCIENTIF / 0428.ZIP / NMR2.BAS < prev    next >
BASIC Source File  |  1985-04-19  |  7KB  |  176 lines
  1. 1 'Program NMR2
  2. 10 'Part 2 of NMRCALC package.
  3. 15 'Calculates and stores eigenvalues and eigenvectors; uses Givens-Householder
  4. 16 ' method of matrix diagonalization.
  5. 20 DEFINT I-N: DEFDBL A-H,O-Z
  6. 30 'COMMON IPFLAG,IREAD,FF$
  7. 31 OPEN "scratch.nmr" FOR INPUT AS #1
  8. 32 INPUT #1, IPFLAG: INPUT #1, IREAD: LINE INPUT #1, FF$
  9. 33 CLOSE #1
  10. 40 DIM A(35,35),V(35,35),E(35),W(5,35),BC(7),FZ(8)
  11. 45 ON ERROR GOTO 60000
  12. 50 CLS:PRINT:PRINT"Now ready to diagonalize the sub-blocks of the Hamiltonian and to get the":PRINT" eigenvalues and eigenvectors.":PRINT: GOSUB 63999
  13. 60 CLS:PRINT:PRINT"Now reading in matrix coefficients to set up reading of sub-blocks.":PRINT
  14. 65 FACTOR = (SQR(5#) - 1#)/2#
  15. 70 DF$ = FF$ + ".inf"
  16. 80  OPEN DF$ FOR INPUT AS 1
  17. 90 INPUT #1,NS
  18. 100 INPUT #1,NF
  19. 110 FOR I = 0 TO NS: INPUT #1,BC(I): NEXT
  20. 120 FOR I = 1 TO NS + 1: INPUT #1, FZ(I): NEXT
  21. 130 CLOSE 1
  22. 140 FOR NQ = 2 TO NS
  23. 150 CLS: PRINT: PRINT "Processing sub-block for Fz =";FZ(NQ)
  24. 160 DF$ = FF$ + "." + RIGHT$(STR$(NQ),LEN(STR$(NQ))-1)
  25. 170 OPEN DF$ FOR INPUT AS 1
  26. 180 INPUT #1, N
  27. 190 FOR I = 1 TO N
  28. 200 FOR J = 1 TO N
  29. 210 A(I,J) = 0
  30. 220 NEXT
  31. 230 NEXT
  32. 240 FOR I = 1 TO N
  33. 250 INPUT #1, A(I,I)
  34. 260 NEXT
  35. 270 FOR I = 1 TO N-1
  36. 280 FOR J = I + 1 TO N
  37. 290 INPUT #1, A(I,J)
  38. 300 A(J,I) = A(I,J)
  39. 310 NEXT
  40. 320 NEXT
  41. 330 CLOSE 1
  42. 340 GOSUB 50000
  43. 350 PRINT:PRINT "Storing eigenvalues and eigenvectors for sub-block.":PRINT
  44. 360 DF$ = FF$ + "." + RIGHT$(STR$(NQ),LEN(STR$(NQ))-1)
  45. 370 OPEN DF$ FOR OUTPUT AS 1
  46. 380 PRINT #1, N
  47. 390 FOR I = 1 TO N
  48. 400 PRINT #1, E(I)
  49. 410 NEXT
  50. 420 FOR J = 1 TO N
  51. 430 FOR I = 1 TO N
  52. 440 PRINT #1, V(I,J)
  53. 450 NEXT
  54. 460 NEXT
  55. 470 CLOSE 1
  56. 480 NEXT
  57. 490 PRINT:PRINT"Storage of eigenvalues and eigenvectors completed.":PRINT:           GOSUB 63999
  58. 500 OPEN "scratch.nmr" FOR OUTPUT AS #1
  59. 510 PRINT #1, IPFLAG: PRINT #1, IREAD: PRINT #1, FF$
  60. 520 CLOSE 1
  61. 600 CHAIN "NMR3"
  62. 50000 EPS = .0000000001#: NMAX = N: NV = N: IORD = 1: M = N
  63. 50010 NM1 = N - 1: IF N=2 THEN 50220
  64. 50020 'Householder transformation.
  65. 50030 PRINT: NM2 = N - 2: PRINT"Tridiagonalizing the matrix.": PRINT
  66. 50040 FOR K = 1 TO NM2: KP1 = K + 1: W(2,K) = A(K,K)
  67. 50050 PRINT"Step:";K: SUM =0
  68. 50060 FOR J = KP1 TO N: W(2,J) = A(K,J): SUM = W(2,J)^2 + SUM: NEXT
  69. 50070 S = SQR(SUM): IF W(2,KP1)<0 THEN S = -S
  70. 50080 W(1,K) = -S: W(2,KP1) = W(2,KP1) + S: A(K,KP1) = W(2,KP1): H = W(2,KP1)*S:      IF H = 0 THEN 50210
  71. 50090 SUMM = 0
  72. 50100 FOR I = KP1 TO N: SUM = 0
  73. 50110 FOR J = KP1 TO I: SUM = A(J,I)*W(2,J) + SUM: NEXT
  74. 50120 IF I >= N THEN 50150
  75. 50130 IP1 = I + 1
  76. 50140 FOR J = IP1 TO N: SUM = A(I,J)*W(2,J) + SUM: NEXT
  77. 50150 W(1,I) = SUM/H
  78. 50160 SUMM = W(1,I)*W(2,I) + SUMM: NEXT
  79. 50170 U = SUMM*.5/H
  80. 50180 FOR J = KP1 TO N: W(1,J) = W(2,J)*U - W(1,J)
  81. 50190 FOR I = KP1 TO J: A(I,J) = W(1,I)*W(2,J) + W(1,J)*W(2,I) + A(I,J)
  82. 50200 NEXT: NEXT
  83. 50210 A(K,K) = H: NEXT
  84. 50220 W(2,NM1) = A(NM1,NM1): W(2,N) = A(N,N): W(1,NM1) = A(NM1,N): W(1,N) = 0:        GERSCH = ABS(W(2,1)) + ABS(W(1,1))
  85. 50230 FOR I = 1 TO NM1: GG = ABS(W(2,I+1)) + ABS(W(1,I)) + ABS(W(1,I+1)):             IF GG > GERSCH THEN GERSCH = GG
  86. 50240 NEXT
  87. 50250 DEL = EPS*GERSCH
  88. 50260 FOR I = 1 TO N: W(3,I) = W(1,I): E(I) = W(2,I): V(I,NV) = E(I): NEXT
  89. 50270 IF DEL = 0 THEN 50530
  90. 50280 'QR method with origin shift.
  91. 50290 PRINT:PRINT"Solving for eigenvalues.":PRINT
  92. 50300 K=N
  93. 50310 L=K: PRINT"E(";K;") = ";E(K)
  94. 50320 IF ABS(W(3,L-1)) < DEL THEN 50340
  95. 50330 L = L - 1: IF L > 1 THEN 50320
  96. 50340 IF L=K THEN 50410
  97. 50350 WW = (E(K-1) + E(K))/2: R = E(K) - WW: ONE = 1: IF WW<0 THEN ONE = -1
  98. 50360 Z = WW - ONE*SQR(W(3,K-1)^2 + R*R)
  99. 50370 EE = E(L) - Z: E(L) = EE: FF = W(3,L): R = SQR(EE*EE + FF*FF): J = L:           GOTO 50390
  100. 50380 R = SQR(E(J)^2 + W(3,J)^2): W(3,J-1) = S*R: EE = E(J)*C: FF = W(3,J)*C
  101. 50390 C = E(J)/R: S = W(3,J)/R: WW = E(J+1) - Z: E(J) = (FF*C + WW*S)*S + EE          + Z: E(J+1) = C*WW - S*FF: J = J + 1: IF J<K THEN 50380
  102. 50400 W(3,K-1) = E(K)*S: E(K) = E(K)*C + Z: GOTO 50310
  103. 50410 K = K-1: IF K>1 THEN 50310
  104. 50420 'Straight selection sort of eigenvalues.
  105. 50430  IF IORD < 0 THEN SORTER = -1 ELSE SORTER = 1
  106. 50440 J = N
  107. 50450 L=1: II=1: LL=1
  108. 50460 FOR I = 2 TO J: IF (E(I) - E(L))*SORTER > 0 THEN 50480
  109. 50470 L=I: GOTO 50490
  110. 50480 II = I: LL = L
  111. 50490 NEXT
  112. 50500 IF II=LL THEN 50520
  113. 50510 WW = E(LL): E(LL) = E(II): E(II) = WW
  114. 50520 J = II-1: IF J>=2 THEN 50450
  115. 50530 PRINT: IF M=0 THEN RETURN
  116. 50540 'Inverse-iteration for eigenvectors.
  117. 50550 QFN = N: EPS1 = SQR(QFN)*EPS: SEPS = SQR(EPS):                                  EPS2 = (.000001#*GERSCH)/(QFN*SEPS): RN = 0: RA = EPS*FACTOR: IG = 1
  118. 50560 FOR I = 1 TO M: PRINT"Processing eigenvector";I
  119. 50570 IM1 = I - 1
  120. 50580 FOR J = 1 TO N: W(3,J) = 0: W(4,J) = W(1,J): W(5,J) = V(J,M) - E(I): RN = RN+ RA: IF RN >= EPS THEN RN = RN - EPS
  121. 50590 V(J,I) = RN: NEXT
  122. 50600 FOR J = 1 TO NM1: IF ABS(W(5,J)) >= ABS(W(1,J)) THEN 50630
  123. 50610 W(2,J) = -W(5,J)/W(1,J): W(5,J) = W(1,J): T = W(5,J+1): W(5,J+1) = W(4,J): W(4,J) = T: W(3,J) = W(4,J+1): IF W(3,J)=0 THEN W(3,J)=DEL
  124. 50620 W(4,J+1) = 0: GOTO 50650
  125. 50630 IF  W(5,J)=0 THEN W(5,J) = DEL
  126. 50640 W(2,J) = -W(1,J)/W(5,J)
  127. 50650 W(4,J+1) = W(3,J)*W(2,J) + W(4,J+1)
  128. 50660 W(5,J+1) = W(4,J)*W(2,J) + W(5,J+1): NEXT
  129. 50670 IF W(5,N)=0 THEN W(5,N)=DEL
  130. 50680 ITERE = 1
  131. 50690 IF ITERE = 1 THEN 50730
  132. 50700 FOR J = 1 TO NM1: IF W(3,J)=0 THEN 50720
  133. 50710 T = V(J,I): V(J,I) = V(J+1,I): V(J+1,I) = T
  134. 50720 V(J+1,I) = V(J,I)*W(2,J) + V(J+1,I): NEXT
  135. 50730 V(N,I) = V(N,I)/W(5,N): V(NM1,I) = (V(NM1,I) - V(N,I)*W(4,NM1))/W(5,NM1):       VN = ABS(V(N,I)): VNI = ABS(V(NM1,I)): IF VNI > VN THEN VN = VNI
  136. 50740 IF N=2 THEN 50790
  137. 50750 K = NM2
  138. 50760 V(K,I) = (V(K,I) - V(K+1,I)*W(4,K) - V(K+2,I)*W(3,K))/W(5,K)
  139. 50770 VNI = ABS(V(K,I)): IF VNI>VN THEN VN=VNI
  140. 50780 K = K-1: IF K>=1 THEN 50760
  141. 50790 S = EPS1/VN
  142. 50800 FOR J = 1 TO N: V(J,I) = V(J,I)*S: NEXT
  143. 50810 IF ITERE > 1 AND VN > 1 THEN 50840
  144. 50820 ITERE = ITERE + 1: IF ITERE > 5 THEN 50830 ELSE GOTO 50700
  145. 50830 'Transformation of eigenvectors.
  146. 50840 IF N=2 THEN 50920
  147. 50850 FOR J = 1 TO NM2
  148. 50860 K = N-J-1: IF A(K,K)=0 THEN 50910
  149. 50870 KP1 = K + 1: SUM =0
  150. 50880 FOR KK = KP1 TO N: SUM = A(K,KK)*V(KK,I) + SUM: NEXT
  151. 50890 S = -SUM/A(K,K)
  152. 50900 FOR KK = KP1 TO N: V(KK,I) = A(K,KK)*S + V(KK,I): NEXT
  153. 50910 NEXT
  154. 50920 J = IG
  155. 50930 IF ABS(E(J) - E(I)) < EPS2 THEN 50960
  156. 50940 J = J+1: IF J<=I THEN 50930
  157. 50950 J = I
  158. 50960 IG = J: IF IG=I THEN 51040
  159. 50970 'Reorthogonalization.
  160. 50980 FOR K = IG TO IM1: SUM =0
  161. 50990 FOR J = 1 TO N: SUM = V(J,K)*V(J,I) + SUM: NEXT
  162. 51000 S = -SUM
  163. 51010 FOR J = 1 TO N: V(J,I) = V(J,K)*S + V(J,I): NEXT
  164. 51020 NEXT
  165. 51030 'Normalization.
  166. 51040 SUM = 0
  167. 51050 FOR J = 1 TO N: SUM = V(J,I)^2 + SUM: NEXT
  168. 51060 SINV = 1/SQR(SUM)
  169. 51070 FOR J = 1 TO N: V(J,I) = V(J,I)*SINV: NEXT
  170. 51080 NEXT
  171. 51090 RETURN
  172. 60000 PRINT: BEEP: PRINT"Error encountered!  Cannot continue.  Will return to main I/O routine.":PRINT"Be sure that required files exist.": GOSUB 63999
  173. 60010 CLOSE 1
  174. 60020 CHAIN "nmr1"
  175. 63999 IF IPFLAG = 1 THEN RETURN ELSE PRINT:INPUT"Hit <Return> to continue.",A$         :PRINT:RETURN
  176.